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© Method for operating a frequency relay and frequency relay for carrying out the method. 

© Method for operating a frequency relay and frequency relay for carrying out the method which, via 
measuring transformers, senses the frequency of an electric power network, determines the value of the 
frequency and delivers an output signal when the frequency thus determined drifts out of a permissible 
frequency range. Specific to the invention is the method of determining the currently prevailing frequency of the 
power network. The measurement signal (y(t)) obtained from the network is transformed, after filtering (1) and 
digitizing (2), into an analytic model (4) in the form of a truncated Fourier series (y(t)), the coefficients of which 
are determined in a parameter estimator (3) which operates with a method of estimation (5) in accordance with 
the least squares method. Starting from model values, according to the invention a computation of the frequency 
can be carried out in a frequency estimator (6), the output signal of which, on the one hand, is returned as 
current frequency value to the parameter estimator and, on the other hand, constitutes a measure of the current 
frequency («) when determining the limits of permissible frequency variations. 
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Method for operating a frequency relay and frequency relay for carrying out the method 

The invention relates to a method for operating a frequency relay according to the precharactertsing 
part of Claim 1. The invention also relates to a frequency relay for carrying out the method. 

Frequency relays are used in a plurality of protective functions in fields like electric power distribution 
or power transformation. Purely generally, it can be said that the duty of such protective devices is to 
5 initiate actions via frequency relays which prevent any harmful effect of a frequency change (rise or 
reduction) going beyond certain predetermined limit values in relation to the nominal frequency of the 
power network in question. 

If. in the event of faults of various kinds, large power generating units (generators or power plants) are 
tripped, the remaining power generation capacity may become smaller than the actual load. This results in 
jo a frequency drop in the network which may grow to such an extent that other gen erators in the network can 
no ionger function. These will then be tripped by their respective under-frequency protection devices. 

An under-frequency protection device is also used to protect the turbine wheels in steam turbines 
against too Sow a speed. In the event of a fault on a turbine governor, an over-frequency protection device 
can influence the stop valve of the turbine and hence prevent the turbine from running at too high speeds. 
75 A frequency protection device can also be used to monitor the speed of synchronous motors and 
rotating reactive power compensators. 

An important object in connection with frequency relays is to determine the currently prevailing 
frequency of the monitored network. Characteristic of the prior art methods for frequency determination is 
that they are based on measuring the distance between the zero crossings of the alternating measurement 
20 signal. However, these methods do not take into account the fact that the mean value of the measurement 
signal may be different frcm zero, which is true above all in connection with and immediately after large 
abrupt changes of the operating conditions of the network. A consequence of this is that the frequency is 
not determined or estimated as accurately as would be of interest in many applications. 

The technique of determining the frequency of a network with the aid of the zero crossings of the 
25 measurement signal has been realized both using classical analog technique and modern numerical 
technique. 

The analog technique comprises several different embodiments. Since the frequency determination 
according to the invention deals with the application of numerical technique, the classical, analog methods 
will be described very briefly and schematically. 

30 in an analog method, a comparison is made between the phase position of the current in an LC circuit, 
the resonance frequency of which is to correspond to the expected frequency of the network, for example 
50 Hz, and the phase position of the current in a purely resistive circuit. Another principle of measurement 
is based on a quartz-controlled oscillator. The number of pulses from the oscillator is counted during each 
cycle or half-cycle of the mains voltage and is compared with a reference level which is set for the 

35 frequency value. A further analog principle is described in SE-A-333 969 entitled "Means for measuring the 
frequency difference between two alternating voltages". In this case, saturable reactors are used as 
frequency sensing elements. 

For the analog technique to operate satisfactorily, normally several rows of filters must be used, and the 
frecuency is usually determined by forming a mean value over several cycles. In case of abrupt changes of 

40 the frequency, a decaying direct voltage component and noise are normally imparted to the measured 
signal. Since the technique is. in principle, based on undisturbed signals, the frequency determination will 
not be as accurate as is desirable. Otherwise, the use of analog components involves trimming and drift 
problems which are not present in numerical processing. 

Characteristic of the application of numerical technique to the frequency determination is that the 

45 measurement signal is regarded as a stochastic process and that an analytical model thereof is formed. The 
numerical technique permits the use of statistical methods for the signal processing, which results in the 
following advantages: a good noise reduction; a faster good Quality measurement result; freedom from 
aging and calibration; no drift problems. 

In the same way as the analog technique can be carried out in many different ways, there are at 

50 present several different ways of carrying out frequency determination using numerical technique. Char- 
acteristic of these, however, is that - in the same way as the analog technique - they are based on 
measuring the distance between the zero crossings of the measurement signal. In exactly the same way as 
for the analog solutions, especially in case of abrupt changes where the mean value of the signal is different 
from zero, this means a certain inaccuracy in the frequency determination. 

Numerical technique applied in this connection is described in a number of publications and patents. 
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inter alia, in M.S. Sachdev and M.M. Giray, "A Least Error Squares Technique for Determining Power 
Systems Frequency", IEEE Trans, on Power Apparatus and Systems, Vol. PAS -104, No. 2, February 1985. 
pages 437-444. The limitations of the described method originate from the assumption that the signal has 
only a fundamental frequency and that oniy the frequencies which lie in the immediate vicinity of the 
5 fundamental frequency are studied. Around this frequency the model is approximated by using Taylor 
series expansions of the trigonometric expressions. A consequence of this is that the accuracy in the 
frequency estimation will be reduced as the actual frequency becomes gradually different from the 
fundamental frequency. 

Otherwise, Sachdev et al make use of the least squares method, which is based on estimating relevant 
70 parameters with signal values obtained from a finite time window. In this way the pseudo-inverse is formed, 
which is a "heavy" process from the computational point of view. 

Frequency determination using numerical techniques based on measurement of the distance between 
the zero crossings of the signal, but using a different technique from that described in the Sachdev et al 
article, is clear from, inter alia, A. A. Girgis and R.G. Brown, "Application of Kaiman Filtering in Computer 
75 Relaying", IEEE Trans, on Power Apparatus and Systems, Vol. PAS -100, No. 7, July 1981, pages 3387- 
3397, and from A.A. Girgis, "A new Kaiman Filtering Based Digital Distance Relay". IEEE Tran. on Power 
Apparatus and Systems, Vol. -101, No. 9, September 1982, pages 3471-3480. 

Concepts from the signal theory, modelling of dynamic systems and numerical analysis have long been 
used within modern feedback control technique to be able to calculate large and complicated process 
20 control systems. To a certain extent these concepts have also started being used in connection with 
frequency determination in frequency relays. The model in this connection often consists of a truncated 
Fourier series expansion of the monitored signal. In the above-mentioned publications by Girgis as well as 
in US-PS-4,455,612 entitled "Recursive Estimation in Digital Distance Relaying System", recursive Kaiman 
filter technique is applied to estimate the unknown parameters of the model. Similar methods have been 
25 used in the previously mentioned article by Sachdev et al and in an article by A.G. Phadke, J.S. Thorp and 
M.G. Adamiak, "A New Measurement Technique for Tracking Voltage Phasors, Local System Frequency, 
and Rate of Change of Frequency", IEEE Tran. on Power Apparatus and Systems, Vol. PAS -102, No. 5. 
May 1983, pages 1025-1038. 

Especially the latter publications deal with investigations in which adaptive parameter estimation 
30 techniques are applied. These techniques work with estimating the parameters over a sliding, fixed time 
window, as also mentioned above, and with the pseudo-inverse used to solve the linear equation system 
which defines the parameters of the model, for example the Fourier coefficients. 

The publication "Theory and Practice of Recursive identification, by L. Ljung and T. Soderstrom. 
Cambridge, Mass., MIT Press, 1983, pages 323-327, describes a method in which parameter adaptation is 
35 recursively applied at each sampling instant with the aid of an exponential, so-called forgetting factor. The 
object of the method, described there is to weigh together old information with an exponentially diminishing 
measure. In this way slow changes are treated. This publication also investigates algorithms for adaptive 
systems with respect to rapidity, convergence and numerical accuracy. 

An abrupt change of the measurement signal unavoidably leads to major changes of the parameters. 
40 For that reason, methods which take care of rapid changes must also be coupled together with the 
parameter estimation, which is otherwise, pure generally, performed in a so-called parameter estimator, to 
determine when an event has occurred in the studied signal. The methods therefore are also described in 
the Ljung and Soderstrom publication, pages 54-60. 

The prior art comprising the normally used Fourier model, which will be described in greater detail. 
45 assumes that the frequency of the studied signal is known with great accuracy, and when studying the 
mains frequency it is assumed that the frequency is 50 Hz or 60 Hz. If this is not the case, the parameters 
in the parameter estimation algorithm will vary periodically with a frequency which is equal to the difference 
in frequency between the frequency in the model and the actual system frequency, which will therefore also 
be true of the estimation errors (see text below) for the studied signal. The amplitude of the oscillation is 
50 also dependent on the frequency difference. This will, of course, entail problems since the model will then 
oscillate because the frequency of the measurement signal, the variation of which is to be measured, 
deviates from the nominal frequency. 

The model technique and the recursive estimation technique, which will now be described in more 
detail, is a summary of the prior art which has been used individually in various applications. The 
55 description and the mathematical expressions and algorithms utilized can be regarded as a concentrate of 
necessary information obtained, for example, from the Ljung and Soderstrom publication, which is needed 
to realize the invention (see pages 321-369). 

If the signal which is to be processed is harmonic with the frequency fa. i.e. can be expressed as 
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y (t + t 9 > = y (t) d) 

In this case a suitable modei of the signal is a truncated Fourier series. This is applicable since voltages 
and currents in a power line system are normally cyclical in stationary state. Also a possible transient phase 
can then be modelled in the same way with sufficient accuracy by assuming that the mean value of the 
signal curve exhibits an exponentially decreasing shape after a fault. The following relation for the model is 
thus valid: 

N 

y (t) = a 0 exp (-bot) + £T c 3 sin (cj 3 + d 3 ) (2) 

This expression is well-known from many applications and. therefore, a description or the quantities and 
parameters included is not necessary. 

Theoretically, M, i.e. the number of terms in the truncated Fourier series, may assume arbitrarily large 
values. Because of the necessity to keep the amount of computation work below a moderate level during 
real-time implementation, N is normally limited so that only the fundamental tone and the first harmonic 
components are included in the processing. 

The accuracy of the estimated parameters is, of course, dependent on the truncation of the Fourier 
series and hence also on the adaptation thereof to the signal. As will be clear from the further description, it 
is also very important that the signal is allowed to pass through an analog filter before the A/D (analog-to- 
digital) conversion is accomplished. 

Since the mathematical notation in connection with signal theory and recursive technique differs from 
the conventional mathematical notation, a few examples will be given to facilitate the reading of the 
following description. 

If 6 constitutes a column vector. 6 T means the same vector in row form. i.e. the transpose of 8. 

If 3 means a true measured value of a vector, then ? means an estimated value of the same vector. 

By "arg" prefixed to a vector function is generally meant an arbitrary vector in the definition set of the 
function. By "arg min" is thus meant that vector which gives the vector function a minimum. 

According to current technology, the model equation for the curve shape of the signal according to 
equation (2) can be transformed to 
y (t) = 6 T o{X) (3) 

where S T represents a vector in row form for estimation of the parameters included in equation (2). i.e. 

e T = (a 0 , -c 0 b 0 , Ci cos di , Ci sin dt . ... c N cos d*. C N sin d*) (4) 

and <a(t) consists of a regressive vector 

o(t) - (1. t. sin <*> c t, cos w 0 t, ... sin N 0 « 0 t, cos N 0 u> 0 t) (5) 

The parameter vector e thus contains elements which provide information about the system whereas <j>- 
(t) comprises the regressors which generate the output signal from the mode!. 

The equations (4) and (5) describe a general model in the sense that they do not require any special 
knowledge of the system save for the periodic structure. As will be clear from equation (5), the elements of 
the regressive vector must be calculated new at each sampling instant. Such a structure can be used in the 
application of the estimation algorithm which will be described below. From (4) it is also clear that the 
following relationships apply between the elements in the 9-vector and amplitude (c,) and phase (dj) in the 
different Fourier components 

c, = (6 2H + eh\) 1,2 (6) 

d, = arg (0 2 2j .„ e 2j ) (7) 

The determination of the system parameters - the parameter estimation - from successively measured 
values of y(t) according to equation (3) is usually carried out according to the least squares method. This 
means that the parameters are chosen to minimize the value of a so-called "loss function" V N . V N is defined 
as follows: 

t 

Vm (A,t) = 3~^- 3 c 3 (j) = *V*,(A,t:-l) + c 2 (t) (8) 
3 = 1 

The function includes a so-called forgetting factor \. which in this connection is a constant, and where «(t) is 
an estimation error function, i.e. eft) = y(t) - y(t). 

The general solution to this equation is given by the estimated value ?(t) = arg mm {V N (\.t)} (9) 
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By differentiating V N with respect tc the parameter vector, which is implicitly included in e(t), the 
following equation (10) is obtained for calculating e(t) which can be solved as soon as the number of 
sampling points exceeds the number of estimated parameters: 



©(t) « /a- - 1 (0) + SZa^-vtj )..--*( j)7 * TL >*-V( j)y( j)l 
L j-l J j-i J 



(10) 



P(0) is equal to the initial value of the covariance matrix of the parameter vector (see latter part of the 
description). For each new updating of the estimated value 5(t), a number of samples are required which 
are equal to or greater than the number of elements in the 0-vector (N > dim 9). This also means that at 
least N values must be sampled before the first estimation can be made unless p~' (0) * 0. 

The estimation of the parameters according to equation (10) can be performed recursively with the aid 
of the foiiowmg algorithm which comprises: 
R(t) = X R(M) + o(t) o r (t) (11) 
R(0) =5'l (12) 
R(t)L(t> = ©it) (13) 
y(t) = 9' (M)o(t) (14) 
e(t) = y(t) - /(t) i'S) 
?(t) = 6 (t-l) - L(t) c(t) (16) 

In addition iz the previously defined quantities, R(t) here designates the covariance matrix of the regressive 
vector,5 its initial parameter, and I a unit matrix. L(t) is a gain vector which is defined according to (13) and 
(19) below. This algorithm contains the solution of a linear equation system which shows that the number of 
arithmetical operations 1 is proportional to N 3 . For real-time implementation, it is therefore more advantageous 
to utilize the low rar,< structure in equation (11) for a recursive real-time algorithm. By applying the matrix 
inversion lemma :o equation (11), an updating of the inverse of R(t) can be obtained, i.e. P(t) = R~'(t), 
whereby, in accordance with known recursive formulae, the following relationships apply: 
r(t) = P(t-1) o ft) (17) 
d(t) = x + o T (:)r(t) (18) 
L(t) = r(t)/d(t) (19) 
P(t) = [P(t-1)-r{t)L 7 m] X (20) 
P(0) = (1-5H (21) 
3(0) = e D (22) 

where d(t) and r(t) denote necessary auxiliary functions. 

The implementation of the algorithm can be carried out in a simple manner since it is exclusively based 
on matrix and vector calculations. The parameter 5 is chosen sufficiently small to avoid problems with 
singular matrices and to enable a rapid adaptation immediately after an initialization. The choice of 0 o can 
either be made arbitrarily, for example all components equal to zero, or empirical model values can be used 
initially until the updating of the covariance matrix has become stabilized. 

For high dimensions of the regressive vector, numerical inaccuracy may cause serious problems. Also 
a limited word length or number of elements tend to have an undesirable effect on the solution; it may, for 
example, become unstable or assume infinite values. Because of these facts it is usually not possible to 
carry out a direct implementation of the equations (17) ... (22). A more appropriate method is then to base 
the solution on an algorithm which is based on square root decomposition of P(t). Such a method is 
described by G.J. Bierman in ''Factorization Methods for Discrete Sequential Estimation", Academic Press, 
New York (1976), pages 44-47. Since P(t) is updated by a low rank structure from P(M), it is possible to 
update the U(t) and D(t) components described in Bierman's publication in accordance with a special 
algorithm mentioned there which reduces the number of numerical operations. This method makes it 
possible to attain the same time complexity as for the conventional algorithm according to equations (17) ... 
(22). For this invention and application, the UD-based algorithm has been investigated with respect to 
numerical accuracy and has proven to be exponentially convergent, which shows that the effects of faults 
introduced into the algorithm, for example as a result of implementation with a curtailed word length, tend to 
55 decay in proportion to the forgetting factor. The point of convergence is, of course, dependent on the 
current number representation. 

The invention aims at developing a method for operating a frequency relay which enables a highly 
accurate, constantly updated and fast determination of the frequency of the input signal which is also 
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accurate under abrupt transient conditions in the monitored power network- 
To achieve this aim the invention suggests a method for operating a frequency relay according to the 

introductory pan of Claim 1. which is characterized by the features of the characterizing part of Claim 1. 
Further developments of the method according to the invention are characterized by the features of the 
5 additional Claims 2 and 3. 

A frequency relay for carrying out the method according to the invention is characterized by the 

features of Claim 4. 

The external functions of a frequency relay according to the invention are the same as those according 
to the state of the art. Specific to the relay, however, is the manner in which the determination of the 

w frequency in question is carried out, namely, using digital technique and microprocessors. This means that 
the frequency measurement as such can be integrated with and be used in other types of modern 
protective relays which require accurate knowledge of the current frequencies of the network. This applies, 
for example, to fault locators and distance protection means of various kinds, the function of which is 
dependent on microprocessors included »n the protection device, and the programming, structural build-up. 

15 etc.. of such microprocessors. 

Immediately after an initialization of the parameter estimation algorithm , the recursion of the covanance 
matrix is sensitive to noise. This effect becomes still more pronounced for small values of X. As an 
immediate consequence of this property, it is possible to find a remedy by applying a time-dependent 
forgetting factor, which is initially equal to a lower value. X Q , and which thereafter converges towards a 

20 predetermined value x<» i.e. 

xft + D = rx(t) + n-r)x» (23) 

X(0) = X 0 (24) 

\(oo) = \co (25) 

where f is a measure of -the rate of convergence towards \«. As usual, suitable choices of x 0 and x» are 
25 given by x o *[0.9.1.0] and fe[0.9.1 .0], respectively. With this tecnnique. a fast convergence as well as good 

numerical properties in cases involving sensitive signals are ootained. 

Otherwise, since in accordance with the prior art, data from a limited time window have been utilized, an 

equation system must be solved at each sampling instant or at least when a new estimation of the 

frequency is desired. This is eliminated by algorithms (17) ... (22). 
30 The invention will now be described in greater detail with reference to the accompanying drawings 

showing - by way of example - in 

Figure 1 a fundamental structure of the frequency determination in a frequency relay according to the 

invention, 

Figure 2 a somewhat more detailed fundamental structure of the frequency determination in a 
35 frequency reiay according to the invention, 

Figure 3 a flow chart for the frequency esiimator included in the invention. 
Figure 4 the sequence of calculations in the parameter estimator included in the invention, 
Figure 5 a structure of a frequency relay according to the invention. 
The principle of the frequency determination according to the invention is clear from Figure 1. To a 
40 certain extent, as also described above, it comprises aspects of the prior art applied in other contexts and 
aspects of new and unique technique. 

The measurement signal y(t) passes through an analog filter 1 and is sampled or converted to digital 
representation by an analog-to-digital converter 2. The digitized signal is passed to a parameter estimator 3 
which has at its disposal a model 4 comprising a structure for forming a model of the measurement signal. 
45 and an estimation unit 5 comprising the structure for the estimation methodology included therein. 

As will have been clear, several different models are available, for example a Taylor series expansion of 
the measurement signal according to Sachdev et ai, or a Fourier series expansion, possibly with an 
exponential mean vaiue according to (2). The model used in this connection is based on a Fourier 
expansion of the signal. 

50 As will aJso have been clear from the above description, several models are available for estimating the 
system parameters or the coefficients in the applied model. The Taylor series expansion of the measure- 
ment signal according to Sachdev et af makes use of a least squares method for determining the 
parameters. In the Girgis publications mentioned above, a technique is used which is based on Kalman 
filtering. 

55 In an application according to the invention, the combination of Fourier series expansion and recursive 

least squares method is used to determine the coefficients of the Fourier series. 

This makes it possible, with the aid of the parameter estimator, the model and the estimation 
methodology, to obtain an estimated, filtered measurement signal y(t) expressed in accordance with the 
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analytical expression of the model. 

As will have been clear from the above, according to current technology an assumption is made as 
regards the frequency and if the frequency of the measurement signal deviates from this assumption, the 
system parameters and hence the residual function e(t) will oscillate with the differential frequency and also 

5 the amplitudes of the residual functions tend to oscillate with the same frequency. If the frequency could be 
estimated and fed into the parameter estimator and the model, this oscillating drawback could be 
eliminated. A condition for this, however, is that the frequency estimation is accurate. 

The method which is comprised by the invention permits an exact calculation of the frequency of the 
measurement signal based on sampled vaiues of the estimated measurement signal, i.e. of y(t). This 

70 calculation, which is accomplished in a frequency estimator 6, is continuously updated. By feeding back the 
estimated frequency to the parameter estimator and the modei, in accordance with elementary feedback 
control technique, the parameter estimator will, after a first frequency assumption, work with the current 
frequency. In consequence thereof, also the output signal of the frequency estimator constitutes an 
estimated value which with great accuracy shows the current frequency of the measured signal. 

75 The invention entails that the Fourier components and the exponential term are uodated recursively 
starting from the immediately preceding estimation. This means that the computational work is considerably 
reduced, which therefore also means that the sampling time can be chosen shorter. The system is built up 
with a fixed number of sampling instants per cycle. The distance between the nodal points is changed as 
soon as a major change of the frequency has been observed. 

20 Since the computation as such can be regarded as a digital filter, no additional software filtering is 
needed. 

The method according to the invention permits the mean value of the measurement signal to be 
estimated after the occurrence of a fault. A consequence of this is that exponentially decaying transients will 
be under control and that reliable estimations of the frequency can be obtained more rapidly. The number 
25 of estimations per cycle is S-2 if S denotes the number of sampling instants during one cycle. The intensity 
of estimations wiil thus be considerably higher with maintained accuracy. 

The need to have knowledge of transients and stationary frequencies is also present in other fields 
where protective relays are used, for example in distance protection devices. Since knowledge of both of 
the above properties is available in the method according to the invention, this knowledge can be 
30 advantageously applied directly to distance protection devices as well. 

In the following a calculation of the frequency of the filtered and modeled measurement signal is 
described: 

If it is assumed that the processed measurement signal can be modelled in accordance with equation 
(2), the stationarily periodic fundamental component of the signal can be extracted by studying 
35 y e (t) = a ' sinwt (26) 

This model contains three unknown quantities, namely amplitude a, angular frequency o> and running 
time t. By sampling three consecutive values y: , y 2 , y3 from the model, which are separated by the 
distance h in time, it is possible to estimate a and v according to the following technique. From equation 
(26) is is clear that the sampled values can be expressed as yi = y e (t-h) = a° (sin(a>t)*cos(cjh) - cos(<jt> 
40 'sin(a>h)) (27a) 
y 2 = y e (t) (27b) 

y3 = y e (t + h) - a°(sin(^t) # ccs(a)h) + cos(«t)'sin(u>h)) (27c) 

By introducing b = a 2 and applying straightforward trigonometric relationships, yt and ya can be 
developed further by substituting y 2 as follows 
45 y. = y 2 cos(u>h) + (b-yi ) t/2 sin(«h) (28a) 
y 3 = y 2 cos(o>h) - (b-yi ) ,/2 sin(u;h) (28b) 
or in a more concentrated form 
yT = -a • sin(u>h-<*>) (29a) 
y 3 - a 4 sin(uh-<*>) (29b) 
50 <p = arc sin(y 2 /a) (29c) 

If equation (28) is solved with respect to sin(ojh), the following is obtained 
sin(a>h) = b-'tyatb-y? ) l;2 -yt(b-yl ) 1 *' 2 ) (30a) 
sin(«h) = b- , (y 3 (b-y| ) 1;2 - y 2 (b-y§ ) 1 ' 2 (30b) 

By introducing x = sin(uh) and manipulating equations (30a) and (30b), the following is obtained 
55 (yi+y 3 )(b-yl )" 2 = y 2 {{b-yi ) vz + (b-y? >" 2 ) (30c) 

The equation only contains one unknown variable, namely b. By continued proceeding, the square roots can 
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be eliminated and the final expression, which is quadratic in b, is given as 
b <«(y« +y3) 2 - 2yl ) 2 -4y 2 )b + 4y \ (yiy 2 ((ys + y') - 2y 1 My 3 + y 
Equation (3D has two solutions. One is zero and the other is given by 



+ y ? ) y 2 )) = 0 



(31) 



4yf (y*Y3-yf ) 



(yx+y3) 2 -4y 2 2 




yxy 3 -yf 



yi-2y 2 +y 3 



y* 



(32) 



The last equality in (32) is favourable if the signal can assume iarge values. Equation (32) is now to be 
interpreted as an estimation of the square of the amplituoe. 

By sertmg z = (30a) = (30b) (33) 
and by introducing in this auxiliary function a computed value of b according to (32). the following 
estimation of the unknown angular frequency is obtained 
x - h~" arc sin(z) (34) 

The mitiai parameters to compute the angular frequency are y-, yz, y: and b. Since the filtered 
amplitude of the s«gnai can be obtained from the parameter estimation algorithm, the computation of b in 
(32) can be eliminated. This means that only the frequency computation according to (34) need be added to 
the original estimation orocess. 

The describee teenmcue utilizes information from three sampled values only. This means that stringent 
demands must be imccsed on the analog filter which precedes the estimator. A better method should 
therefore be to consicer 1 . a sequence of measured values <y- . j = 1:N) and to successively compute 
estimations cf the angu«a: frequency in accordance with 
xn = bp 1 fyn(b.-y£ ) } 2 - Yn.» (bn-y n >"" 2 ) (35) 

It is assumed nere :~.a: b„ is obtained from the parameter estimator. It is clear from equation (35) that 
only one square root neec be computed if |b n -bj<i? for some small number t>. 

To prevent too great deviations from the true value, a sliding mean value is used as the final estimation, 
Z(t), of the angular frequency 
w(t) = a°u>(t-1) + n-ui • x(t) (36) 

The filter parameter jl is chosen here in the same interval as the previously described forgetting factor 
X. To obtain a rooust estimation, one additional precautionary measure should be taken. If the estimation of 
the frequency computec at a given sampling instant differs too much, for example by p 0 V from the running 
mean value, it shouid be omitted. It is also advantageous to omit such estimations which deviate sufficiently 
from the immediately preceding one. To limit the number of incorrect values, the filtered values are used in 
place of the instantaneous ones. 

The methoc described requires an accurate, analog pre-filtering of the measured signal values in order 
not to render the frequency estimation too inaccurate. In addition, the number of overlooked estimations 
must be kept at a limitea level if the studied signal is regular. 

From equation (28) it is possible to derive the following expression for computation of the angular 
frequency, given three consecutive sample values of the proceeded signal 



J5 y^ys 

h- 1 arc cos (37) 

2y a 



which can also be expressed as 

Z - h" [arc sm(y 3 a) - arc sin(y 2 a)] = h- 1 (arc sin [b-'(y 3 (b-yf )' 2 - y 2 (b-y| ?'*)] (38) 

To calculate u> according to these relationships, however, trigonometric expressions must be developed. 
For this reason, it is aovantageous to derive a method which reduces the number of mathematical functions 
to be evaluated. Such an equation is provided by (35). 

The signal modei usea in equation (26) is limited to the special case where o = 0. However, this is not 
a serious restriction since the phase angle is a relative quantity which is only of significance in cases where 
more than one signal is considered. In other cases the omission of the phase angle in the model should not 
influence the frequency estimation other than during the transient phase. By taking the phase angle into 
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account, the algorithm can be extended to be applicable also to the estimation of the relative phase 
distance between different signals. On the other hand, this information is already available in the Fourier 
model. In order to make the representation complete, also the extended algorithm is presented here. 
Let the signal model be represented by (2) and define y2 as 
5 y e (t) = y(t) -a 0 e-^ (39) 

Contrary to the previous case, the function in equation (36) is developed at four consecutive sampling 
instants ti-h, ti, ti +h and U + 2h to gain knowledge of the phase shift <p. By analogy with (30) the following 
is obtained: 

sin («h + o) = b-'(y 2 (b-yJ? ) i;2 -yi(b-yi ) 1 ' 2 ) (40a 
10 sin (wh + o) = b- , (y 3 (b-y| ) 1/2 - ya(b-y§ ) 1/2 (40b) 
sin i«h + o) * b-'Jyitb-yl ) 1 ' 2 - y3(b-y| ) 1/2 ) (40c) 

From equations (39) and (40) the amplitude and the argument uh + <t> can be estimated with an analog 
technique which has been described above, that is to say, 
sin(w + 0) = c (41a) 
75 sin(2ajh + 0) = sin(ci>h-o)cos(«h) + cos(wh + <*>)sin(wh) = d (41b) 
A combination of the expressions in (41) gives 
c • cos (wh) + (1-c 2 ) 1 ' 2 sinwh = d (42) 

Simple trigonometric relationships make it possible to represent (41) in the more succinct form 
sin (ujh + a) = a (43a) 
20 sin (a) = c (43b) 

This means that the angular frequency can be estimated with the aid of the following relationship u> = 
h~" (arc sin(d) - arc sin(c)) (44) 

and in this way an estimation of the phase angle can be obtained as 

o = arc sin(c) - uh = 2arc sin(c) - arc sin(d) (45) 
25 By using (41) the variables c and d can be expressed as a function of the sampled values according to 

c = b- , (v 2 (b-y? ) 12 - yi(b-yi ) v2 (46a) 

d = b-'(y4(5-y§ ) ,/2 -y3(b-y2 )i* (46b) 

Depending on the system to which the estimator is to be applied, it may be necessary to study further 

expansions of the algorithm (44). After a fault the mean value is initially different from zero and it will 
30 decrease exponentially towards 2ero according to equation (2). Also when an extra signal source has been 

connected, the mean value is different from zero. In the latter case the mean value remains different from 

zero even after the transient phase, which means that it is necessary to estimate it. However, the mean 

value is estimated in the main algorithm, as has been described previously. By taking this into account, the 

frequency estimator can be applied to oscillating signals which are superimposed on the mean value. From 
35 this fact the conclusion may be drawn that it is advantageous to apply one frequency estimator for each 

phase. 

In the following the technical basis for frequency determining embodiments will be described: 
It is clear from equations (31), (34) and (36) that it is necessary upon each updating of the frequency 
estimation to carry out the following operations: (additions/subtractions, multiplications, divisions, taking of 

40 square roots) = (4,7,1,2). The square roots can be determined either by us ing the Newton-Raphson 
algorithm, or by carrying out a linear interpolation in a pre-computed table. In the former case the maximum 
number of operations will be (a,m.d) = (10,10,10). This corresponds to 5 iterations, which is an acceptable 
number if the initial solution is sufficiently good. To obtain as good an algorithm as possible with a minimum 
number of iterations, it may be advantageous to scale the argument so that it lies in the interval (1,10). This 

45 means that 10 1/2 and 0.1 ,;2 must be pre-computed and saved for subsequent treatment of the argument if it 
has been scaled. The total number of arithmetical operations for an updating, with the exception of scaling, 
will thus be (a,m,d) = (24,27,21). If a TMS32010 is used as processor during the implementation, the time 
of execution will be (24 + 27 + 21 * 4) • 3us = 405txs. This consideration assumes that a division takes 4 
times as long as a multiplication. A few additional operations must be carried out during the scaling, which 

so means that 0.5 ms is a likely time of execution for one estimator. It should also be considered that one 
parameter estimator for each phase is required. On the other hand, it is sufficient with only one frequency 
estimator. The shortest sampling time will therefore not exceed 1 ms. From this follows that 20 samples per 
cycle can be carried out, that is to say, 19 independent estimations can be carried out. This is not too short 
a sampling time for the estimator, since the changes in the fundamental frequency are mostly rather small. 

55 i.e. ± 1 Hz/s. A faster estimation procedure can be obtained by linear interpolation in a square root table 
where the nodal points of the argument are equally spaced. 

A very important matter when implementing the frequency estimator is to reduce the number of 
operations as much as possible. It is aiso necessary to apply such algorithms that are convergent and 

9 



;SDOClD: <EP 0297418A1 J_> 




EP 0 297 418 A1 



robust so that they can be implemented with a short word length. The time of execution for one it eration of 
the algorithm aepends essentially on the following factors: 1) the number of estimated parameters. 2) the 
word length in question, 3) the algorithms in question. 4) the cycle time of the microprocessor, 5) the 
implementation of the algorithm, for example several parallel activities if possible, and 6) rendering the 
5 algorithm discrete during the numerical processing. Also the accuracy is dependent on some of the listed 
points. In addition, the correspondence in appearance between the modei and the real system plays a 
decisive role. 

It is clear from equations (17) ... (22) that if x and u> assume constant values, the gain vector L(t) will 
stationariiy only depend on these parameters. This means that the L(t) vector can be pre-computed as a 

70 first step before the parameter estimation algorithm is initialized. From this follows that the number of 
operations is reduced from being quadratically dependent on N. i.e. the number of terms in the truncated 
model, to being linearly dependent on N. On the other hand, the need for data memory is at the same time 
increased. If X and <d are limited to a finite number of discrete values, the data memory can be kept at a low 
level. In actual fact, the forgetting factor can normally be maintained constant during the whole execution 

?5 instead of prescribing a time dependence or values which are pointwise constant. Since there are rigorous 
demands to keep the frequency within a small interval around the nominal frequency, only a limited amount 
of ^-values will be of interest. If the set of possible values lies sufficiently close to each other, it is possible 
to use linear interpolation between the nodal points of w. A special processor can be connected to the 
system tc measure the frequency and to calculate a new set of gain vectors if changes of the prescribed 

20 magnitude have occurred. If it is suitable, the last step can also be a part of the initialization procedure. 

In this connection it is also of interest to investigate after how many cycles L(t) can be considered to be 
stationary. It can be shown that the deviation is less than 1CT 2 after 5 cycles for x = 0.98 and 50 samplings 
per cycle, i.e. altogether 750 sampling instants. For X = 0.95, which corresponds to 3 cycles, in total 450 
sampling instants must be saved for each frequency. If the computer memory is not sufficient to save the 

25 gain vector for a high value of X . it will be necessary to successively upaate the covariance matrix P(t) 
according to (10) and to compute L(t) from (20). However, the amount of computations increases 
consideraoiy and for such an implementation, of course, a lower sampling frequency must be chosen. 

When the regressive vector is formed, cos(nw 0 t) and sin(nw 0 t) must be evaluated. As in the case of the 
gain vector, this can be performed either via linear interpolation in a table or by using a fast algorithm for 

30 trigonometric functions. The first alternative is. of course, the fastest and should thus be chosen if sufficient 
data memory is available to store the tables. 

Another important property of the regressive vector (5) proves to be that the updating of the covariance 
matrix is independent of the signal to which the estimator is applied. This fact makes it possible to limit the 
handling to one gain vector. When applying the algorithm to a three-phase system, the parameters in the 

35 different phases can then be updated according to 
Sx(t) = e x (t-1) + L(t) [x(t) - e l (t-1) o(t)] (47) 
where x = l A , Ib. lc> U A , U B , U c . respectively. 

In the following preferred embodiments of the invention wilt be described: 
Figure 2 shows roughly how, in principle, the building blocks in the frequency determining part of a 

40 frequency relay are connected, purely from the functional point of view. The figure is thus a somewhat more 
detailed specification of the blocks 3,4.5 and 6 from Figure 1 . 

The incoming, filtered and sampled measurement signal y(t) is supplied to a building block or module 7 
comprising an identification algorithm, and to a module 4 comprising the model in question. The 
identification algorithm which comprises the equations (17) ... (22) is used to produce numerical values of 

45 the parameters in the given model according to equation (2) and the processed signal. As will have been 
clear from the foregoing description, the combination of the given model and the identification algorithm 
assumes that the fundamental frequency is known with high accuracy. It is therefore necessary, according 
to the invention, to connect an estimator 8 for estimation of the frequency. The frequency estimator, which 
is based on the equations (34), (35) and (36), is supplied with the numerical values from the identification 

so algorithm, and by cooperation with the modei, the model frequency is successively adjusted to the 
variations in the frequency of the studied signal. 

The regressive vector according to equation (5) is. in principle, generated in module 9 at each sampling 
instant to constitute an input quantity to the identification algorithm. 

The event handling according to module 10 comprises the study of the transient process after a rapid 

55 and sufficiently great change in the system dynamics with the aid of the loss function defined according to 
equation (8). 

As soon as the value of the loss function exceeds a certain value, the parameter estimation algorithm is 
re-initialized via module 1 1 and equations (21 ) and (22). 
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Fed-out quantities via module 12 consist of the estimation vector d and the estimated fundamental 
frequency Z or ?. 

Figure 3 illustrates roughly a flow chart of how the algorithm for the frequency estimation can be 
implemented. All the quantities included are assumed to be represented as a fixed point number with 15 
5 bits and one bit for the sign. When implementing the algorithm, due account is taken of the fact that the 
variables are of different orders of magnitude by choosing scale factors during the arithmetical operations. 

The task and function of the different blocks are clear from Figure 3 as well as from the following 
description of the blocks: 

Block 13 comprises a starting function, 
/o Block 14 comprises 

(i) initialization of the following parameters: 
X forgetting factor 
n f number of harmonic components 
n s number of samples per cycle 
75 n p number of cycles to be stored 
ft last approved estimated frequency 
tg sampling time 

p pointer to pre-computed gain vector 
5 initialization parameter for covariance matrix 
20 KC number of samples to be used to start up covariance matrix and stabilization of the parameter 
vector, standard setting = 20 

NF number of samples to be moved after re-start of parameter estimator before frequency estimator is 
switched in, standard setting = 20 

(jj) generation of a table of trigonometric functions over a quarter cycle. 
25 (iii) initialization of the parameter vector d. If no advance information of the system exists, al! 

elements are set equal to zero. 

Block 15 comprises initialization of the covariance matrix as P(0) = l/S, where I, as mentioned before, 
is a unit vector. 

Block 16 comprises resetting of the initialization counter. 
30 Block 17 comprises adjusting the initialization counter upwards at each sampling during start-up of 

the covariance matrix. 

Block 18 comprises parameter estimation according to Figure 4. 
Block 19 ensure that the start-up of parameter estimator proceeds for NC samples. 
Block 20 comprises resetting of the frequency updating counter. 
35 Block 21 comprises adjusting the frequency counter upwards. 

Block 22 comprises storage of the last sampled value of the signal as y(t), as is clear from the figure. 
Block 23 comprises, in the same way as block 18, parameter estimation according to Figure 4. 
Block 24 comprises filtering of the signal by the estimated parameters according to (14). A window of 
the last three filtered signal values is stored in the memory to be used during the frequency estimation. 
40 Block 25 checks whether the time t from the last initialization of the parameter estimation is great in 

relation to the time constant r of the exponential term. 

Block 26 reduces the number of estimated parameters when the exponential term has decayed. 
Block 27 ensures that the frequency is held constant during the first NF updatings after an 
initialization. 

45 Block 28 checks whether V(\,t) has a large amplitude. If V exceeds a certain value ci. an abrupt 

change in the signal has occurred, which justifies carrying out an immediate restart of the parameter 
estimator. 

Block 29 comprises producing an estimated value of the frequency according to (34). 
Block 30: If the estimated frequency f e at a given sampling instant deviates by more than « 2 from the 
so last approved estimate of the frequency, the latter is updated by the former. Otherwise, the next sample is 
carried out immediately. 

Block 31 stores the computed frequency if it deviates to a sufficient extent from the immediately 
preceding accepted frequency estimation. 

Block 32 computes a new sampling time t s - (n s • f..)"' so that the number of samplings per cycle is 
55 constant. 

Block 33 generates a tripping signal as output quantity form the protective relay. 
A more detailed description of the parameter estimation according to blocks 18 and 23 is clear from 
Figure 4 wherein 
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block 34, i.e. "START-UP", means that the parameter estimator is to be initialized, 
block 35 means that the next sample value form the A/D convener is read. 

block 36 comprises the formation of the regressive vector. Linear interpolation between the nodal 
points is used when the argument deviates too much from the nodal points. Trigonometric relationships are 
5 utilized to calculate regressors which correspond to harmonics, 

block 37 comprises calculation of the gain vector L(t) according to equation (20), 
block 38 updates the covariance matrix Pit) according to equation (19). 

block 39 updates the estimation vector 9(n) according to equations (15) and (16) and updates the loss 
function V(\.t) as equation (8), and 

10 biock 40 stops the parameter estimator. 

The embodiments described so far have substantially comprised the functional part of the frequency 
measurement part of a frequency relay according to the invention. Figure 5 illustrates an electrical/physical 
structure of the frequency relay. To the monitored network RST there is connected either a current 
transformer. CT, or a capacitive voltage transformer CVT. designated as 41. As in Figure i, Figure 5 

75 comprises an anaiog filter 1 and an analog-to-digital converter 2. The A/D converter 2 digitizes the analog 
signal upon a call from a first microprocessor I. 42. In this connection, it is sufficient with a resolution of 12 
bits. 

The described algorithms for frequency and parameter estimation can be implemented with the aid of 
the above-mentioned first microprocessor 42 and a second microprocessor II. 43. One of these may 

20 advantageously be a signal processor. The processors are working towards a common RAM/ROM memory 
44. The execution of the processors is performed simultaneously. In processor I the parameters estimator is 
implemented in accordance with equations (17) ... (22) and (10). Processor II is used to estimate the 
frequency according to equation (34) while at the same time the updating of the parameters takes place in 
processor l. This means that the execution in processor II is time-displaced by one sampling interval 

25 compared with processor I. Also logic handling and initialization of the algorithms are carried out in 
processor II. 

The implementation of the various algorithms can, of course, be carried out in a many different ways 
oepenoing on the available hardware and the invention is therefore not restricted to the described realization 
as shown m :he figures. 

30 The estimated frequency u> or f obtained is then compared in a conventional manner in a comparator 45 

with the lower permissible limiting frequency f- and the corresponding upper frequency f u and a TRIP signal 
is obtained when any of the limits is exceeded. 

If the frequency estimation is to be carried out in a three-phase system, the availability of the three 
measurement signals should be utilized so that the frequency estimation works towards the first voltage in 

35 the order R.S.T. which is not zero. This can be performed by studying the Fourier coefficients for the 
voltages which are obtained with parameter estimators for the different phases. Only one generation for the 
regressive vector is necessary since it can be used for all the parameter estimators. The computation of the 
gain vector is carried out in a separate processor. One processor for each phase is used for the parameter 
estimation and the logic handling. This means that tow estimation vectors are updated in each processor. A 

40 supervisory processor is used for administration. 

Claims 

45 1 . Method for operating a frequency relay which is used for protecting objects in a power network which 

are sensitive to frequency variations, said frequency relay receiving as an input signal, via current or voltage 
transformers connected to the power network an analog measurement signal y(t) containing information 
about the frequencies occurring in the power network, said frequency relay further being designed to deliver 
an output signal when the frequency (f) of the network drifts out of a given frequency range, the method 

50 comprising a sampling and digitizing of the measurement signal, a model generation of the measurement 
signal in the form of a truncated Fourier series and an exponential part as well as a recursive least squares 
estimation method, characterized in that the recursive least squares estimation methodology is used for 
estimation of the parameters in the model of the measurement signal, that the frequency in question is 
computed with the aid of a number of model values consecutively filtered and sampled from the model, and 

55 that the value of the frequency thus computed is consecutively supplied to the recursive least squares 
estimation methodology and is used as a base for the parameter estimation and the frequency computation. 
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2. Method according tc Claim 1, characterized in that the frequency o> of the model is computed with 
the aid of three sampled, filtered and consecutive model values yi,y 2 and y 3 , which are separated by a 
time interval h in accordance with the following equation 

^ = h-* arc cos (37) 

2y 3 



3. Method according to Claim 1, characterized in that the frequency w of the model is computed witn 
10 the aid of four sampled, filtered and consecutive model values y1.y2.y3 and y*. which are separated by a 
time interval h according to the following equation 
cj = h- 1 • (arc sin(d) - arc sin(c)) (44) 
and that the phase angle <t> is estimated as 
0 = arc sin(c) - wh = 2arc sin(c) - arc sin(d) (45) 



T5 



wherein 



C = 



20 



25 



b = 



b-*[<y4(b-y§ - y3(fc-y*) 1/2 ] 

YxY3-yf 



4y 2 
yi+2y a +y3 



yx-2y 2 +y3 



(46a) 
(46b) 

(32) 



4. Frequency relay for carrying out the method according to any of the preceding claims, which 
30 frequency relay, via current or voltage transformers (41) connected to the power network, is adapted to 
receive as input signal an analog measurement signal (y(t)) containing information about the currently 
prevailing frequency of the power network, the frequency relay further being designed to deliver an output 
signal (TRIP) when the frequency (f) of the network is found to lie outside a given frequency range, 
characterized in that the frequency relay comprises an analog filtering device (1). an analog-to digital 
35 converter (2), microprocessors (42,43) for data processing in connection with a model generation of the 
measurement signal, estimation of the parameters of the model and estimation of the currently prevailing 
frequency based on modei values and other event handling of the data processing, a RAM/ROM memory 
(44} and a comparison device (45) adapted to evaluate whether the computed frequency (f) lies within 
permissible limits, i.e. is greater than a lower limiting frequency (fi) or smaller than an upper limiting 
40 frequency (f u ), and to deliver a TRIP signal when these limits are exceeded. 



45 



50 



55 
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